Background

This code contains model files and associated material for the 2023 ICES course on Stock Synthesis held in Copenhagen, Denmark.

Model setup

# install.packages("devtools")
# devtools::install_github("r4ss/r4ss", ref="development")
# 
# install.packages("caTools")
# library("caTools")
# 
# install.packages("r4ss")
library(r4ss)
library(here)
#remotes::install_github("PIFSCstockassessments/ss3diags")
library(ss3diags)
library(knitr)
library(tidyverse)
dir1<-here("herring-aspm")# version 3.30.21
dir2<-here("Workshops" , "SS3_ICES_Course", "herring-just-2019-as-equilibrium") 
dir3<-here("herring-scaa") # 
dir4<-here("Workshops" , "SS3_ICES_Course", "herring-scaa-2dar") # 
dir5<-here("herring-scaa-2dar-noTV") 
dir6<-here("Workshops" , "SS3_ICES_Course", "predators") # v. 3.30.21

Run models

#shell(cmd="ss_win") # run SS para windows
# aqui debo cambiar el setwd() para correr el modelo
#dir5 <- setwd("~/DOCAS/SA_Krill/s5")
# r4ss::run(dirvec = dir3, model = './ss_osx', 
#                     skipfinished = FALSE) 
# #system('./ss_osx')

#another way to run with new r4ss

r4ss::run(
  dir = dir1,
  exe = "ss_osx",
  skipfinished = FALSE, # TRUE will skip running if Report.sso present
  show_in_console = TRUE # change to true to watch the output go past
)

r4ss::run(
  dir = dir1,
  exe = "ss_osx",
  skipfinished = FALSE, # TRUE will skip running if Report.sso present
  show_in_console = TRUE # change to true to watch the output go past
)

r4ss::run(
  dir = dir2,
  exe = "ss_osx",
  skipfinished = FALSE, # TRUE will skip running if Report.sso present
  show_in_console = TRUE # change to true to watch the output go past
)

r4ss::run(
  dir = dir3,
  exe = "ss_osx",
  skipfinished = FALSE, # TRUE will skip running if Report.sso present
  show_in_console = TRUE # change to true to watch the output go past
)

r4ss::run(
  dir = dir4,
  exe = "ss_osx",
  skipfinished = FALSE, # TRUE will skip running if Report.sso present
  show_in_console = TRUE # change to true to watch the output go past
)

r4ss::run(
  dir = dir5,
  exe = "ss_osx",
  skipfinished = FALSE, # TRUE will skip running if Report.sso present
  show_in_console = TRUE # change to true to watch the output go past
)

r4ss::run(
  dir = dir6,
  exe = "ss_osx",
  skipfinished = FALSE, # TRUE will skip running if Report.sso present
  show_in_console = TRUE # change to true to watch the output go past
)

Diagnosis

herring-aspm Model

Saco los outputs en html

Data disponible para este escenario. Espinel es la serie mas consistente del conjunto de datos.

SSplotData(base.model1, subplot = 2)
## Warning: The `subplot` argument of `SSplotData()` is deprecated as of r4ss 1.45.1.
## ℹ Please use the `subplots` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Respecto a los valores y parametros biologicos modelados, los siguientes graficos identifican los estimadores puntuales del recurso

SSplotBiology(base.model, 
              subplots =1, 
              labels = c("Length (cm)", "Age (yr)",
                         "Maturity", 
                         "Mean weight (kg) in last year",
                         "Spawning output",
                         "Length (cm, beginning of the year)", 
                         "Natural mortality",
                         "Female weight (kg)", 
                         "Female length (cm)", 
                         "Fecundity", 
                         "Default fecundity label",
                         "Year", 
                         "Hermaphroditism transition rate",
                         "Fraction females by age at equilibrium"),
 )

aporte de las cohortes por año para las capturas.

SSplotCohortCatch(base.model1)

AJuste de tallas por flota

SSplotComps(base.model1, subplots = 1)

Otros plots

SSplotDynamicB0(base.model1, uncertainty = F)
## Plotting Dynamic B0

#SSplotSPR(base.model3)
SSplotPars(base.model1)
## Plotting distributions for 4 estimated parameters (deviations not included).
## $mcmc not found in input 'replist', changing input to 'showpost=FALSE'

Salida de las biomasas con las dos flotas

SSplotTimeseries(base.model1, subplot = 1)

SSplotTimeseries(base.model1,
  subplot=1,
  add = FALSE,
  forecastplot = FALSE,
  uncertainty = FALSE,
  bioscale = 1,
  minyr = -Inf,
  maxyr = Inf,
  plot = TRUE,
  print = FALSE,
  plotdir = "default",
  verbose = TRUE,
  btarg = "default",
  minbthresh = "default",
  xlab = "Year",
  labels = NULL,
  pwidth = 6.5,
  pheight = 5,
  punits = "in",
  res = 300,
  ptsize = 10,
  cex.main = 1,
  mainTitle = FALSE,
  mar = NULL
)

herring-scaa Model

Saco los outputs en html

Data disponible para este escenario. Espinel es la serie mas consistente del conjunto de datos.

SSplotData(base.model3, subplot = 2)

SSplotRecdevs(base.model3)

Respecto a los valores y parametros biologicos modelados, los siguientes graficos identifican los estimadores puntuales del recurso

SSplotBiology(base.model3, subplots =1, 
              labels = c("Length (cm)", 
                         "Age (yr)", 
                         "Maturity", 
                         "Mean weight (kg) in last year",
                         "Spawning output",
                         "Length (cm, beginning of the year)",
                         "Natural mortality",
                         "Female weight (kg)",
                         "Female length (cm)", 
                         "Fecundity", 
                         "Default fecundity label",
                         "Year", 
                         "Hermaphroditism transition rate", 
                         "Fraction females by age at equilibrium"),
               )

aporte de las cohortes por año para las capturas.

SSplotCohortCatch(base.model3)

AJuste de tallas por flota

SSplotComps(base.model3, subplots = 1)

Otros plots

SSplotDynamicB0(base.model3, uncertainty = F)
## Plotting Dynamic B0

#SSplotSPR(base.model3)
SSplotPars(base.model3)
## Excluding 57 deviation parameters because input 'showdev' = FALSE
## Plotting distributions for 14 estimated parameters (deviations not included).
## $mcmc not found in input 'replist', changing input to 'showpost=FALSE'

Salida de las biomasas con las dos flotas

SSplotTimeseries(base.model3, subplot = 10)
SSplotTimeseries(base.model3,
  subplot=1,
  add = FALSE,
  forecastplot = FALSE,
  uncertainty = FALSE,
  bioscale = 1,
  minyr = -Inf,
  maxyr = Inf,
  plot = TRUE,
  print = FALSE,
  plotdir = "default",
  verbose = TRUE,
  btarg = "default",
  minbthresh = "default",
  xlab = "Year",
  labels = NULL,
  pwidth = 6.5,
  pheight = 5,
  punits = "in",
  res = 300,
  ptsize = 10,
  cex.main = 1,
  mainTitle = FALSE,
  mar = NULL
)

herring-scaa-2dar-noTV Model

Saco los outputs en html

Data disponible para este escenario. Espinel es la serie mas consistente del conjunto de datos.

SSplotData(base.model5, subplot = 2)

Respecto a los valores y parametros biologicos modelados, los siguientes graficos identifican los estimadores puntuales del recurso

SSplotBiology(base.model5, subplots =1, 
              labels = c("Length (cm)", 
                         "Age (yr)", 
                         "Maturity", 
                         "Mean weight (kg) in last year",
                         "Spawning output",
                         "Length (cm, beginning of the year)",
                         "Natural mortality",
                         "Female weight (kg)",
                         "Female length (cm)", 
                         "Fecundity", 
                         "Default fecundity label",
                         "Year", 
                         "Hermaphroditism transition rate", 
                         "Fraction females by age at equilibrium"),
               )
## Note: this model uses the empirical weight-at-age input.
##       Plots of many quantities related to growth are skipped.
## Warning in max(biology[["Len_mean"]]): no non-missing arguments to max;
## returning -Inf

aporte de las cohortes por año para las capturas.

SSplotCohortCatch(base.model5)

AJuste de tallas por flota

SSplotComps(base.model5, subplots = 1)

Otros plots

SSplotDynamicB0(base.model5, uncertainty = F)
## Plotting Dynamic B0

#SSplotSPR(base.model3)
SSplotPars(base.model5)
## Excluding 69 deviation parameters because input 'showdev' = FALSE
## Plotting distributions for 3 estimated parameters (deviations not included).
## $mcmc not found in input 'replist', changing input to 'showpost=FALSE'

Salida de las biomasas con las dos flotas

SSplotTimeseries(base.model5, subplot = 10)
SSplotTimeseries(base.model5,
  subplot=1,
  add = FALSE,
  forecastplot = FALSE,
  uncertainty = FALSE,
  bioscale = 1,
  minyr = -Inf,
  maxyr = Inf,
  plot = TRUE,
  print = FALSE,
  plotdir = "default",
  verbose = TRUE,
  btarg = "default",
  minbthresh = "default",
  xlab = "Year",
  labels = NULL,
  pwidth = 6.5,
  pheight = 5,
  punits = "in",
  res = 300,
  ptsize = 10,
  cex.main = 1,
  mainTitle = FALSE,
  mar = NULL
)

Desembarques

Los desembarques utilizados para cada una de las flotas que generan remoción en el recurso reineta, a saber; espinel, enmalle e industrial (Figura 6).

\label{F2}Desembarques de merluza común por flota

Desembarques de merluza común por flota

Predator fleet with RW

Random Walk (RW) refers to a mathematical model that describes a stochastic process in which a variable changes randomly over time, without a clear trend or pattern.

Specifically, a random walk can be used as a Bayesian estimation technique to infer the posterior distribution of an unknown parameter. In this approach, it is assumed that the prior distribution of the parameter is a normal distribution with mean zero and a known variance, and that the parameter value at each time step follows a random walk process. Based on the observed data and the prior distribution, the posterior distribution of the parameter can be calculated using Bayesian inference.

Random walk is a useful tool for parameter estimation in dynamic models, as it allows for modeling uncertainty and variability in the parameter’s evolution over time. However, it is important to note that the random walk assumes that the changes in the parameter are random and without a clear trend, which may not be appropriate in all cases.

#########################################################################
#COMPARACION DE MODELOS con distinto desembarque
#########################################################################


#PLOT labels, name of each model run
legend.labels <- c('herring-aspm',
                   'herring-just-2019-as-equilibrium')

#read in all model runs
#note if cover=T you need a hessian; if covar=F you do not need a hessian
biglist <- SSgetoutput(keyvec = NULL, 
                       dirvec = c(dir1,dir2,dir3),  
                       getcovar = F)

#create summary of model runs from list above
summaryoutput <- SSsummarize(biglist)

SSplotComparisons(summaryoutput, subplots = 1:20, plot = TRUE,
                  print = T, endyrvec = "default", indexfleets = NULL,
                  indexUncertainty = FALSE,
                  indexQlabel = TRUE, 
                  indexQdigits = 4, 
                  indexSEvec = "default",
                  indexPlotEach = FALSE, 
                  labels = c("Year", "Spawning biomass (t)",
                             "Relative spawning biomass", 
                             "Age-0 recruits (1,000s)",
                             "Recruitment deviations", 
                             "Index", 
                             "Log index", 
                             "1 - SPR", 
                             "Density",
                             "Management target", 
                             "Minimum stock size threshold",
                             "Spawning output",
                             "Harvest rate"), col = NULL, 
                  shadecol = NULL, pch = NULL,
                  lty = 1, lwd = 2, spacepoints = 10, staggerpoints = 1,
                  initpoint = 0, tickEndYr = TRUE, shadeForecast = TRUE,
                  xlim = "default", ylimAdj = 1, xaxs = "r", yaxs = "r",
                  type = "o", uncertainty = TRUE, 
                  shadealpha = 0.1, legend = TRUE,
                  legendlabels = "default", legendloc = "topright",
                  legendorder = "default", legendncol = 1, sprtarg = NULL,
                  btarg = NULL, minbthresh = NULL, 
                  pwidth = 6.5, pheight = 5,
                  punits = "in", res = 300, ptsize = 10,
                  plotdir = "~/DOCAS/Workshops/SS3_ICES_Course", 
                  filenameprefix = "", 
                  densitynames = c("SSB_Virgin", "R0"),
                  densityxlabs = "default", densityscalex = 1,
                  densityscaley = 1, densityadjust = 1, 
                  densitysymbols = TRUE,
                  densitytails = TRUE, 
                  densitymiddle = FALSE, 
                  densitylwd = 1,
                  fix0 = TRUE, new = TRUE, 
                  add = FALSE, 
                  par = list(mar = c(5, 4, 1, 1) + 0.1), 
                  verbose = TRUE, mcmcVec = FALSE,
                  show_equilibrium = TRUE)



SStableComparisons(summaryoutput,
                   likenames = c("TOTAL", 
                                 "Survey", 
                                 "Length_comp", 
                                 "Age_comp", 
                                 "priors",
                                 "Size_at_age"),
                   names = c("Recr_Virgin", 
                             "R0", "steep", 
                             "NatM",
                             "L_at_Amax", 
                             "VonBert_K", 
                             "SSB_Virg", 
                             "Bratio_2017", "SPRratio_2016"),
                   digits = NULL, modelnames = "default", csv = FALSE,
                   csvdir ="~/DOCAS/Workshops/SS3_ICES_Course", 
                   csvfile = "parameter_comparison_table.csv", verbose = TRUE,
                   mcmc = FALSE)

Analisis retrospectivo

#do retrospective model runs
SS_doRetro('~/DOCAS/Workshops/SS3_ICES_Course', 
           '', newsubdir = "retrospectives",
           subdirstart = "retro", years = 0:-5,
           overwrite = TRUE, exefile = "ss",
           extras = "-nox -nohess", 
           intern = FALSE, CallType = "system",
           RemoveBlocks = FALSE)

retroModels <- SSgetoutput(dirvec=file.path('~/DOCAS/Workshops/SS3_ICES_Course',
                                            "retrospectives",
                                            paste("retro",0:-5,sep="")))

retroSummary <- SSsummarize(retroModels)
endyrvec <- retroSummary$endyrs + 0:-5

SSplotComparisons(retroSummary, 
                  endyrvec=endyrvec, 
                  legendlabels=paste("Data",0:-5,"years"),
                  plotdir='~/DOCAS/Workshops/SS3_ICES_Course',
                  plot=TRUE,print=T)

TableCompare <- SStableComparisons(retroSummary,
                                   likenames=like, 
                                   names=names, 
                                   modelnames=c('B','-1','-2','-3','-4','-5'), 
                                   csv=TRUE, 
                                   csvfile="RetroRuns.csv",verbose=FALSE)

Course Notes

  • Lorenzen mortality works just in case you have parametric growth. Used a lot in US in Europe (Cardinale, 2023). About F:
  • Both use the same F calculations
  • Hybrid best when F is low
  • At high F, hybrid’s algorithm to maintain match to catch inhibits movement towards convergence
  • Parameter approach can be configured to get starting value for F parameter during early phases using hybrid, then converts to parameter to finish the run
  • F method = 4 gives best flexibility because it allows hybrid vs parameter conversion on a fleet-specific basis (today is the best method. If you have fleet specific).
  • 0.23 # F ballpark value in units of annual_F. Is usefuul when we have outlier data (High biomas or high F) and you can set what year happen this, i. e. ; -2008 # F ballpark year (neg value to disable).
  • Body weigth is in kilograms -Options 3 & 4 are experimental in SS3, but match some other models, like small pelagics stock.
  • If fixing without having an intended rationale, then let the estimated recruitments deviate from the curve: –- use recdev option 2, not 1, for more flexibility –- set lambda on recdevs to low value (<0.10)